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We investigate kinetic pathways of the DNA melting transition using variable-range versions of 
the Poland-Scheraga (PS) and Peyrard-Dauxois-Bishop (PDB) models of DNA. In the PS model, 
we construct a 0'^-field theory to calculate the critical droplet profile, the initial growth modes, and 
the exponent 7 characterizing the divergence of the susceptibility near the spinodal. In the PDB 
model, we use a mean field analysis to calculate 7. We compare these theoretical results with Monte 
Carlo and Brownian dynamic simulations on the PS and PDB models, respectively. We find that 
by increasing the range of interaction, the system can be brought close to a pseudospinodal, and 
that in this region the nucleating droplet is diffuse in contrast to the compact droplets predicted by 
classical nucleation theory. 



PACS numbers: 

I. INTRODUCTION 

The DNA melting transition is an interesting theoret- 
ical problem because a quantitative understanding of its 
mechanism may provide insight into how biological en- 
zymes physically interact with DNA [T]. While much 
theoretical work has examined the role of large nonlinear 
excitations as a precursor to melting [2H1], surprisingly 
little focus has been paid to kinetics of the transition 
itself. In this article, we use classical and spinodal nu- 
cleation theory in conjuncton with simulations to study 
the kinetics of melting in both short- and long-range ver- 
sions of the Peyrard-Dauxois-Bishop (PDB) [3l |4j and 
Poland-Scheraga (PS) \, 5J models of DNA. 

Nucleation plays an important role in many systems 
undergoing a phase transition [6 . In homogeneous sys- 
tems, a droplet forms from a spontaneous fluctuation 
and grows into the stable phase. Before nucleation, the 
system is trapped in a metastable well in a free energy 
landscape. The system samples the phase space of this 
metastable well until a fluctuation drives the system to 
the top of a free energy barrier that separates the sta- 
ble and metastable wells. At this point, the droplet is 



referred to as critical and the system is equally likely 
to nucleate or return to the metastable well. DNA is 
believed to undergo such a process because its sharp 
melting curve indicates a first-order phase transition [T] 
and because hysteresis in the melting curve suggests the 
existence metastable states |7HS]. In addition, possible 
nucleation bubbles have been observed via electron mi- 
croscopy |T0] . 

There is experimental evidence for the presence of long- 
range (LR) interactions in DNA pi lTTHTB] . Telestability 
experiments on block copolymers show cooperativity ef- 
fects over at least 10-15 base pairs (bp) [12] while differ- 
ential melting curves with multistep behavior show co- 
operatively melted regions of 100-350 bps [T]. In low 
salt concentration, cooperatively melted regions can be 
as a large as a few thousand bps [1 1] . Both experimen- 
tal and molecular dynamics studies suggest that bases 
beyond nearest-neighbor can effect the enthalpic change 
that arises from the opening of a given bp [1^1 [IT] • In ad- 
dition, it is known that nearest-neighbor (NN) PS models 
underestimate the probability of single bp opening 119) . 
Only recently have researchers begun examining how LR 
interactions affect the dynamics of DNA models [50] . To 
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consider the effects of the interaction range on the ki- 
netic pathway, we simulate the Peyrard-Dauxois-Bishop 
(PDB) [3 g and PS [H] models with both NN and 
LR interactions. 

It is well known that LR systems undergoing a phase 
transition can be quenched into metastable states near a 
pseudospinodal pTM23j . In the mean field (MF) limit 
i? — oo, this pseudospinodal becomes a well-defined 
spinodal 23J, which is the limit of metastability. The 
barrier to nucleation for a d-dimensional system scales 
with the interaction rang e i? as i?'^A/i3/2-'i/4^ where Ah 
is the distance away from the spinodal 23 . As such, 
the nucleation rate is much smaller for LR systems, and, 
practically, one must get very close to the spinodal to 
observe nucleation. Nucleation near a spinodal or pseu- 
dospinodal is similar to classical nucleation near the coex- 
istence curve in that droplets become critical by reaching 
the top of a free energy barrier and then either grow or 
decay with equal probability. Unlike classical nucleation, 
which is initiated by compact droplets that resemble the 
stable phase [23] , spinodal nucleation is characterized by 
the formation of diffuse fractal-like droplets whose am- 
plitude differs little from the metastable background [251 - 
[27] . In addition, the growth modes of classical droplets lie 
on the droplet surface, while spinodal droplets grow from 
their center. For these reasons, the inclusion or exclusion 
of LR interactions in DNA models makes a substantial 
difference in regard to the character of nucleation that 
can be observed. 

In addition to changing the qualitative shape of the 
critical droplet, moving the system toward a pseudospin- 
odal also causes the isothermal susceptibility x to di- 
verge. This effect has been observed in supercooled wa- 
ter [28| . and its measurement may be a useful method 
for determining whether or not nucleation in real DNA 
exhibits spinodal effects. While there is inherent fuzzi- 
ness in the definition of the pseudospinodal in real sys- 
tems [29] , practically one can determine both its location 
and the exponent 7 characterizing the divergence if the 



metastable lifetime is longer than the measurement time 
of X- If £^ spinodal exists for DNA, it should in priniciple 
be possible to measure both the spinodal temperature Tg 
and 7. 

Finally, it is important to note that biological DNA 
is an intrinsically heterogeneous system because of its 
pseudo-random sequence of bases. In this work, we con- 
centrate on homogeneous nucleation, which is purely ini- 
tiated by spontaneous fiuctuations, rather than in het- 
erogeneous nucleation, where nucleation is aided by the 
presence of boundaries, defects, or other impurities [30j . 
However, since heterogeneous sequences of bases are 
clearly important biologically, we simulate random base 
pair sequences to determine what effects the inhomo- 
geneities have on the kinetic pathways of the transition. 

This paper is organized as follows. Sec. II provides a 
detailed description of the modified PDB and PS models. 
In Sec. Ill, we postulate a phenomenological free energy 
for the PS model and use this to calculate the shape of 
the droplet profile, initial growth modes, and the suscep- 
tibility exponent. In Sec. IV, we provide a MF calcula- 
tion of the susceptibility exponent in the PDB model. In 
Sec. V, we describe the results of our simulations of long 
and short range PDB and PS models including evidence 
for nucleation, the shape of droplet profiles and growth 
modes, and the calculation of the spinodal exponent. We 
summarize and discuss the practical implications our re- 
sults in Sec. VI. 

II. PDB AND PS MODELS 

We modified both the PDB and PS models to include 
long range interactions. In the modified PDB model, the 
Hamiltonian is given by 

N .2 

i=l 

where the state of each bp is specified by its separation 
Hi and its time derivative iji . The first term is the kinetic 
energy of bps with combined mass m — 300 amu and the 
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second term is a potential given by 

y ({y J) = A If + E (2) 

The first term in Eq. [2] is an on-site Morse poten- 
tial describing the net attraction between strands due 
to a combination of hydrogen bonding, solvent interac- 
tions, and the repulsion of negatively charged phosphate 
groups O m EI] . For the homogeneous case the dissocia- 
tion energy Di is treated as a constant D = 0.04 cV while 
for the hctcrogeneoues case it is given by Dat — 0.032 
eV and Dqc = 0.048 cV for A-T and C-G bps, respec- 
tively. The inverse well width is given by a = 4.45 A^^. 
The term W{yn,ym) represents the stacking interaction 
and is given by the anharmonic potential 

where K = 0.06 eV/A, a = 0.35 A-\ and p = 0.5. 
In the original model, this term is purely NN such that 
i? = 1 . We allow this interaction to extend to range 
i? > 1 to examine the effect LR interactions have on the 
kinetic pathways of the transition. 

In the modified PS model, statistical weights are given 
to bound and unbound segments. A bound segment is 
energetically favored because of hydrogen bonding and 
stacking interactions. Unbound segments are entropi- 
cally favored because single stranded loops have a much 
shorter persistence length allowing them to sample a 
larger configuration of phase space. Our LR PS model 
can be described by the Hamiltonian 
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where cr^ = — 1 and Ui — +1 represent bound and open 
bps, respectively. Here, i?o,i represents the binding en- 
ergy which is assumed to be the same for all bps in the ho- 
mogeneous case. In the second term, the parameter Ko 



represents the interaction between adjacent base pairs. 
This term has been added so that the effects of LR in- 
teractions in the model can be studied. In order to have 
roughly consistent parameters between the two models, 
we have set i?o.i = -Di and Kq — K{1 + p). 

The final term in Eq. [4] represents an effective poten- 
tial due to entropic effects caused by to loops of size I. 
Here, we have chosen s — 74.4 to give a biologically rel- 
evant melting temperature T = 350 K. The exponent 
c = 2.15 is chosen consistent with simulation results on 
self-avoiding random walk loops [THIIIS]. Choosing c > 2 
ensures that the melting transition will be first order. 
The cooperativity il = 0.3 is chosen larger than values 
published elsewhere so that small loops may be observed 



III. 



-FIELD THEORY IN THE PS MODEL 



Without the entropic term describing the degeneracy 
of loop configurations, the PS Hamiltonian of Eq. [4] is 
isomorphic to the Ising model in the lattice gas repre- 
sentation. This suggests it may be possible to convert 
the PS model into a (/I'^-field theory in the same way as 
the Ising model. For this reason, we postulate a Landau- 
Ginzburg- Wilson (LGW) free energy functional of the 
form 

= / dr[R^{\/(j)f/2 + K(f>^ + h(f>~ TSi<j))]. (5) 



(6) 



Here, (t>{r) is the coarse-grained magnetization 

1 NT^ 



(r) 



^CG 



and LcG is the size of the coarse-grained region. From 
this definition, it is clear that (/)(r) e [—1,1]. Here, the 
parameter h — Lcg{Eq + Ko)/2 is the energy associated 
with flipping a coarse-grained region against the direction 
of an Ising-like magnetic field. The parameter k sets the 
critical temperature of the field theory. For convenience, 
we chose k = Kq to match the microscopic model. The 
exact numerical choice of parameters is not expected to 
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change the quahtative resuhs. The first three terms in 
the integrand of Eq. [5] describe the energetic contribu- 
tions from the PS model while the last term describes 
the entropic contributions arising both from the coarse- 
graining and from the degeneracy associated with un- 
bound loops. 

In principle, the entropy S{<j)) in Eq. [s] can be calcu- 
lated in a similar way as the entropy that arises when 
coarse-graining an Ising model to obtain a field theory. 
The difficulty of calculating this term a priori stems 
from the entropy's dependence on the distribution of loop 
sizes. To circumvent this problem, we assume the entropy 
S{(f)) of a given coarse-grained region is equal to the aver- 
age entropy of all loop distributions consistent with the 
magnetization (/) of the region. This approximation is ex- 
pected to be reasonable when the loop sizes are smaller 
than LcG- For simplicity, we choose a particular coarse- 
graining size LcG = 16 and enumerate the 2^^ possible 
combinations of spins. For each combination, we cal- 
culate (j) and the loop entropy J2ioops^^S{^s'- /l'^) of the 
configuration. We then calculate the average entropy for 
a given ip, 



\ loops 



(7) 



This numerical determination of S{4>) can be fit to a 
fourth order polynomial 



(8) 



where the coefficients are given by bi = 27.5, 62 = 7.06, 
&3 = 3.83, and = -0.701. 

Combining Eq. [5] and Eq. [8] we now write our LGW 
Hamiltonian as 

"i?2(V0)2 



F{(j)) = / dr 



(9) 



where 

= ^Tb4^^-Tb3cl>^~{Tb2^K)cl>^ + {h~Tb,)^. (10) 

In Fig. [1] we plot the free energy density /{(j)) for 
T = 350 K, T = 410 K, and T = 471.1 K corresponding 




FIG. 1: (color online). Plots of the free energy density / 
vs. (f} for various T. Plots are shown for T — 350 K (black, 
solid), T = 410 K (red, dashed), and T = 471.1 K (blue, 
dash-dotted). 



roughly to the melting temperature T,„, an intermediate 
temperature, and the spinodal temperature T^. In this 
plot, (j) is restricted to lie within a range — 1 < < 1, 
consistent with its definition. At T=:350 K, we see that 
/((/)) has two wells at (/) w 1 and (f) ^ — I corresponding to 
the open and bound states, respectively. Since the melt- 
ing temperature should lie on the coexistence curve, one 
would expect these wells to be the same depth, rather 
than having a deeper bound state well as show in the fig- 
ure. This discrepency arises because of the assumption 
of small loops and will be considered shortly. Near the 
spinodal at T = 471.1 K, the metastable well has almost 
vanished. As suggested in the plot of f{(p) for T = 410 K, 
/((/)) evolves continuously with a gradual disappearance 
of the metastable well as the spinodal temperature is ap- 
proached. 

From the LGW free energy, we calculate the shape 
of the critical droplet. The critical droplet is a saddle 
point in the free energy landscape satisfying the Euler- 
Lagrange equation: 



SF 



Sf 



dr^ 



= 0, 



(11) 



There are two solutions to Eq. [TT] that are independent 
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FIG. 2: (color online), (a) Critical droplet profiles are shown 
for r = 350 K (black, solid), T = 410 K (red, dashed), and 
T — 471.1 K (blue, dash-dotted). A green line has been 
added to show the unphysical region where > 1. (b) The 
initial growth mode v{r) for each temperature of part (a). A 
vertical green line has been added to denote the region where 
the assumption of small droplets clearly fails when T is small. 

of r. One of these represents the metastable bound state 
= a-iid the other, which would correspond to the 
stable unbound state, appears at > 1 and is unphys- 
ical. The actual stable unbound state occurs at cj) ~ 1. 
A spatially nonconstant solution (f>{r) with the bound- 
ary condition that 0(oo) ~ (j)MS represents a fluctuation 
away from the metastable well. This fluctutation is the 
critical droplet. 

Eq. [TT] is analogous to the equation of motion for a 
particle in a potential V{(j)) — ~f{(t>) with r representing 
time |:24 . In this analog, the boundary condition ^ — 
at r = corresponds to the particle having no initial 
velocity. The particle starts up the side of the larger 
hill (stable minimum) and rolls off until finally coming 



to rest on top of the smaller hill (metastable minimum.) 
This equation of motion was solved numerically using a 
fourth-order Runge-Kutta method for T = 350 K, 410 
K, and 471.1 K. The resulting nucleating droplet profiles 
are shown in Fig. [2^. 

A green line has been added to the figure to emphasize 
where (f) > 1. At T ~ 350 K, the peak of the predicted 
critical droplet lies above this line. This unphysical re- 
sult is due to the assumption that loops are much smaller 
than the coarse-graining size Lea = 16. Our method of 
calculating S{(f>) undercounts the entropic contribution 
of loops with size I > 16, because each coarse-grained re- 
gion's entropy is calculated individually without regard 
to loops that may extend into the next coarse-grained 
cell. This discrepency is significant near the coexistence 
curve where large compact droplets are expected to nu- 
cleate the system. Undercounting large loops increases 
the predicted value of the free energy for the unbound 
state, which can be obeserved in the plot of /(0) in Fig. 
[T] At the melting temperature, the free energy of un- 
bound state should be identical to that of the bound 
state, but the theory predicts a significantly higher free 
energy for the dissociated state. Since the magnetization 
is greater than one for a finite region around the cen- 
ter of the droplet, the actual critical droplet should be 
compact in this region, and the free energy of this actual 
compact droplet will be lower than the free energy pre- 
dicted by the field theory which undercounts the entropic 
contribution. 

Away from the coexistence curve, the field theory is 
a more accurate description because nucleating droplets 
are no longer large and compact. From the shape of 
the droplet near at r=471.1 K (Fig. [2] inset), this is ap- 
pears to be the case. If one notes the change in scale for 
the droplet near the spinodal, it is clear that the droplet 
must be diffuse since its amplitude differs little from the 
metastable background. At T = 410 K, (Fig. (2)3), the 
droplet is intermediate between spinodal and classical, 
consistent with previous results [131 [55] . 
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The critical droplet is perched atop a saddle point in 
the free energy landscape. During the initial growth away 
from this point, the system rolls off the saddle along the 
path of steepest descent. If we write fluctuations to the 
critical droplet as (t>{r) = 4>{r) + v{r), then in the neigh- 
borhood of the saddle, the LGW free energy can written 
as F{(j)) = F{4>) + F"{v), where 



F"{v) 



dr 



A^2 



(12) 



The normal modes of the system near the critical droplet 
configuration are solutions to the Schrodinger equation 



(13) 



There is one negative eigenvalue for this equation cor- 
responding to an instability. This initial growth mode 
increases exponentially with time. 

For each of the temperatures listed above, the growth 
modes were calculated numerically using the shooting 
method and are depicted in Fig. [2]3. As before, the field 
theory works poorly near the coexistence curve where 
droplets are large and compact. A green line has been 
added to denote where the T = 350 K droplet crosses 
into the unphysical region in which > 1. In this re- 
gion, the T=350 K droplet is expected to be compact. 
From Fig.[2]D, one can see that the theory predicts maxi- 
mum growth for a r=350 K droplet occurs at the center 
of the droplet. A compact droplet cannot grow from its 
center because this region is already in the stable phase, 
hence v{r) must be zero for small r. Since v{r) ^ for 
small r, the maximum of the growth mode must lie on 
the surface of the droplet when the system is near the 
coexistence curve. Away from the coexistence curve, the 
theoretical results are more accurate since the droplets 
are diffuse. From the figure, one can see the predicted 
position of maximum growth remains at the center of the 
droplet as the spinodal is approached. Unlike droplets 
that form near the coexistence curve, droplets that form 
near the spinodal are diffuse and can grow from their 
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FIG. 3: Log-log plots of x{T) vs. Ts - T for the PS model 
obtained from (a) the field theory and (b) simulations. A 
linear plot of x vs. T is shown in the inset of both figures. 



center. This agrees with previous results on nucleation 
in Ising models |25j . 

Since the spinodal is a critical-like point, the sus- 
ceptibility X — ^^ni^ is expected to diverge as the 
spinodal temperature Tg is approached. We calculated 
X ~ <P{h+Ah)-4)(h) various temperatures within the 
metastable region using Ah — 10^^. As the T —i' Tg, a 
power law divergence of the form x ^ \Ts — Tl^'^ is ob- 
served as shown in Fig. [3^. Fitting this plot to a power 
law, we obtain a T5 « 471 K and 7 « 0.5. 

The value of 7 can be explained if one compares it to 
the analogous divergence x = l^s — in the Ising 

model. The coefficient of the linear term in a cj)'^ map- 
ping of the Ising free energy is the field h, while the lin- 
ear coefficient of the PS free energy is linearly dependent 
on temperature. As such, the temperature divergence 
in the PS model should be the same as the field diver- 
gence in the Ising model, leading to an exponent 7 = 0.5. 



Such must be the case as the T divergence in the Ising 
model characterized by 7 = 1 only occurs by lowering T, 
whereas the spinodal in the PS model is approached by 
raising T. 

IV. MEAN FIELD THEORY IN THE PDB 
MODEL 

While constructing a (p'^-Held theory would be appre- 
ciably more difficult in the PDB case, we can still ob- 
tain reasonable results for the exponent 7 with a simpler 
mean field theory. We start from the long-range Peyrard- 
Dauxois-Bishop model of DNA. Since all base pairs in 
the mean field model are equivalent, we first write the 
Hamiltonian of a single base pair, 

Hsingie {v, {Vj}) = Vm(2/) + Vi {y, {yj}) . (14) 



Here, 



VMiv) = D (e-'^y - ly 



(15) 



is the on-site Morse potential and 

vi{y,{yj}) = §iJ2iy^-^yy^ + y^ 
j 

+pe-"ye-''y^{y^-2yyj+y])} 

is the interaction potential. The sum is over all base 
pairs yj within range R. The kinetic energy term has 
been dropped for simplicity. 

In mean field, we take the limits N ^ 00 and R 
N/2. In the Hamiltonian above, there are five types of 
interaction terms that depend on the separation of neigh- 
boring base pairs. Within the sum, these terms are pro- 
portional to yj, yj, e~°'y^ , yje~"yi , and yje"""^ . In order 
to have a self-consistent mean field theory, it is necessary 
not only that 



but also 



(y') = ^T. 



and 



(y^-"') = ^12yi^~"''' 



N 

J 

All five of these relations must hold for self- consistency. 

To construct a self-consistent mean field theory, we 
first rewrite the mean field Hamiltonian as 

+ K{y'^ - 2yci +C2+ pe-°'y{y'^C3 - 2yci + C5)} 

At present, ci, C2, C3, C4, and C5 will be treated as pa- 
rameters that can take arbitrary values. In order for the 
theory to be self-consistent, these parameters must be 
chosen such that 

ci = (y), 
C2 = {v') , 
C3 = (e-"«), 

C4 = (ye-"^). 



and 



C5 = {y^e-^y) . 



(16) 



In order to compute the self-consistent values for the 
above parameters, we use the following procedure. We 
begin by assigning the values Ci = 0, C2 = 0, C3 = 1, 
C4 = 0, and C5 = 0. Using these starting values, we de- 
termine a new set of parameters from the Maxwell Boltz- 
mann probability distribution P{y) oc exp{— /3i7MF(2/)}- 
This distribution can be normalized by dividing by the 
partition function 

/oo 
dye^^{-pHMF{y)]. (17) 
-00 

Once normalized, this distribution can be used to calcu- 
late new values for the parameters from the relation, 
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(O) = j dyOP{y). 



(18) 
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FIG. 4: (color online). Plots of the free energy F vs. Ci for 
T = 7547^ (black, solid), 1100 K (red, dashed), and 1450 K 
(blue, dotted). For T < 1100 K there is a well near ci — 1. 
At T ~ 1100 K, the well disappears. While the presence of a 
well suggests metastability, this well only appears because we 
have restricted the system to change only along the ci axis. 



Using equation (18), new values for each of the five 



paramerters were computed numerically. This procedure 
was repeated until all five parameters had converged on 
some final value. 

We are particularly interested in finding the spinodal 
temperature. At this temperature the metastable well 
disappears. To find this temperature, we use a simi- 
lar procedure to the one described above. This time, 
we solve self-consistently for the parameters C2, C3, C4, 
and C5 but choose the value of ci. Physically, this corre- 
sponds to allowing all the base pairs to reach equilibirum 
with the constraint that they must maintain a certain 
mean separation. This will only be self-consistent for 
certain values of ci . In Fig. |4] we plot the free energy 
F(ci) = —kBTln Z{ci) as a function of ci for several 
temperatures. Here, Z(c\) is a restricted partition func- 
tion, which is obtained by integrating over states with a 
particular (not necessarily self-consistent) value of ci. 



Z(ci)= / dyP{y). (19) 

J —QfO 

For T < 1100 K one can see what appears to be a 



FIG. 5; (color online). Plots of (y — ci) vs. Ci for T = 
711 K (black, solid), 827 K (red, dashed), and 943 K (blue, 
dotted). For T < 827 K there are two zeros that correspond 
to metastable and unstable fixed points. At T = 827 K the 
two fixed points collide and give a spinodal. 



metastable well. While this figure is useful for illustration 
purposes, it gives the false impression that the spinodal 
occurs at T sa 1100 K. This is incorrect. While there is 
clearly a well when the free energy is projected on the 
Ci axis, we have not considered what happens to the free 
energy when we change any of the other four parameters. 
Unless the well is a local minimum in the space of all five 
parameters, then the system will not be trapped at this 
value of (y). As it turns out, this projection of the free 
energy overestimates the actual spinodal temperature. 

To find the correct spinodal temperature, we again use 
a similar procedure to the one described above. As be- 
fore, we solve self-consistently for the parameters C2, C3, 
C4, and C5 but leave c\ as a free parameter. We then 
calculate the mean separation (y) as a function of c\ . In 
Fig. [5] we plot (y — ci) vs. c\ for several temperatures. 
This plot was obtained by numerically integrating the 
Maxwell Boltzmann distribution, 



(20) 



Self-consistency requires that {y — ci) = 0. Positive val- 
ues of (y — ci) indicate that the mean base pair separa- 
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FIG. 6: Log-log plots of x{T) vs. Tg - T for the PDB model 
obtained from (a) the mean field theory and (b) simulations. 
A linear plot of x vs. T is shown in the inset of both figures. 



tion would be larger than ci in equilibrium as individual 
base pairs vifould drift toward larger separations. Neg- 
ative values of {y — Ci) indicate that the mean separa- 
tion should be smaller than ci because base pairs would 
drift toward smaller values. The first zero in the plot 
is then a stable fixed point corresponding to a stable or 
metastable minimum. When a second zero appears, it 
represents an unstable fixed point. At T w 827 K the 
two zeros fuse. Above this temperature, there are no ze- 
ros, so the system is unstable. As such, this temperature 
corresponds to the limit of metastability and represents 
the spinodal. If one attempts to find {y) for T > 827 
using the full self-consistent, s/he finds that the parame- 
ters will not converge and the mean separation will grow 
without bound. 

Using the position of the fixed point and the defini- 



tion x{T) — ^^f^, we calculated x(T) and the exponent 
7 associated with susceptibility's divergence. We approx- 
imate the derivative as 



d{y) 



dh 



{y)h=dh - iy)h=o 



dh 



(21) 



with dh = 1 X 10^12. From the xiT) vs. T in the inset 
of Fig. [3^, we see a sharp divergence as we raise the 
temperature. Using a two parameter fit, we observe 7 — 
0.5 as r — Ts, consistent with the results for the (f)'^-tie\d 
theory of the PS model. 

V. SIMULATIONS RESULTS 

We simulate the PDB and PS models using Brownian 
dynamics (BD) and the Metropolis Monte Carlo (MC) 
algorithm, respectively. In BD, the system evolves via a 
Langevin equation 



miji = ~WV{yi) - 7?) + 77(0 



(22) 



where the noise ri{t) is random Gaussian with {ri{t)) — 
and {r]{t)r]{t')) = 2jkBTS{t - t'). In time units of 
T = 1.018 X 10~^'*s, we use a time step of 0.25 r and 
a damping constant 7 = 10"** . In the MC simula- 
tions, random spins are flipped and the change in energy 
between the original and final states is calculated. Nega- 
tive changes in energy are always accepted, while positive 
changes are accepted with probability exp(— AE'/Zc^T), 
where /S.E is the change in energy. Both simulations use 
periodic boundary conditions for simplicity. In the next 
two subsections, we describe nucleation in homogeneous 
systems while in the third subsection we consider nucle- 
ation in heterogeneous systems. 

A. Evidence for Metastability 

We plot the time evolution of both models for ranges 
i? = 1 and i? = 205 in Fig. [7] For the PDB model, 
we monitor the growth of the mean separation of the 
strands (y) = ^ yi/N, while for the PS model we plot the 
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FIG. 7: Metastability and Nucleation. Plots are shown for 
the time evolution of (a) the TV = 128 PS model with R = 1 
quenched to T = 365 K, (b) the N = 128 PDB model with 
R = 1 quenched to T = 380 K, (c) the iV = 4096 PS model 
with R = 205 quenched to T = 470 K, and (d) the iV = 4096 
PDB model with R = 205 quenched to T = 800 K. 

"magnetization" (cr) = ^di/iV. Each run is quenched 
instantaneously from T — 300 K to a higher tempera- 
ture listed in the caption. Each temperature is chosen to 
give a metastable lifetime that is long enough to clearly 
demonstrate metastability but short enough to have a 
reasonably fast run time. 

For both ranges in the PS model, we observe a pe- 
riod where the magnetization M stabilizes before grow- 
ing sharply. Growth is visibly sharper for systems with 
long-range interactions. Within the time allotted for the 
runs, the system never returned to the bonded state af- 
ter complete separation. Stability prior to spontaneous 
growth is the hallmark of nucleation out of a metastable 
state. In addition to these signatures, we observe the 
formation of a droplet (see below) that grows into the 
stable phase much like those observed in real DNA. 

In the PDB model, we observe similar stability and 
growth of the mean bp separation (y) for both ranges, 
again with sharper growth for the long-range system. 



Unlike the PS model, the NN PDB model will fre- 
quently rebind even after the strands have completely 
separated. Bases in real DNA bubbles sample a large 
three-dimensional phase space, while bases in PDB bub- 
bles can only move in one dimension. For this reason, 
bases in PDB bubbles are more likely to find and re- 
bind with their complementary pair than bases real DNA 
bubbles which must search a larger space. The result is 
analogous to a one-dimensional random walker which will 
inevitably return to its starting position after some finite 
time. Recombination of the strands was not observed in 
runs of the LR PDB model. While it is possible that 
one might observe recombination of strands in a LR sys- 
tem given enough time, a more likely explanation is that 
the LR interactions create a greater entropic barrier that 
suppresses the likelihood of reforming the double strand. 

It should be noted that both the PS and PDB models 
exhibit noticibily smaller fluctuations when the interac- 
tions are long-range. This is expected since increasing the 
interaction range takes the system closer to its mean field 
approximation in which there are no fluctuations. Prac- 
tically speaking, reducing the size of fluctuations is ben- 
eficial since one need not wait very long to observe that 
the system has reached metastable equilibrium. This is 
particularly useful in long-range systems, where the sim- 
ulation speed is inherently slower. In contrast, the NN 
PDB model exhibits fluctutations large enough that one 
must observe long runs in order to clearly see that the 
system has reached metastable equilibrium. 

Since nucleation is an activated process governed by a 
constant rate, it should follow Poisson statisics. We de- 
termine the nucleation rate for both models and ranges 
by measuring the nucleation time of an ensemble of sys- 
tems run with a different random noise. We define the 
nucleation time as the time at which either the mean 
magnetization in the PS model or the mean separation 
in the PDB model is greater than some threshold. The 
threshold values were chosen to be sufficiently larger than 
metastable fluctuations and are listed in the caption of 
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FIG. 8: Histogram of nucleation times for (a) the PS model 
with iV = 128 and i? = 1 quenched to T = 370 K, (b) the 
PDB model with N = 128 and = 1 quenched to T = 400 
K, (c) the PS model with TV = 4096 and R = 205 quenched 
to T = 470 K, and (d) the PDB model with iV = 1024 and 
R = 50 quenched to T = 780 K. 

Figs. |8j in which we plot on a log-scale histograms of the 
nucleation times. In both long- and short-range models, 
we observe exponential decay for large times, indicating 
that systems reach the stable state at a constant rate. 
Strand separation at a constant rate is consistent with 
the notion that this DNA melting occurs via nucleation. 

B. Divergence of the Susceptibility 

We calculated the susceptibility at various tempera- 
tures by measuring the fluctuations in the PS model, 



XT 

and the PDB model, 

XT = 



T 



T ■ 



(23) 



(24) 



Here, and {■■■)^ signify time and spatial averages, 

respectively. 

As one approaches the spinodal, the system nucleates 
quickly. To suppress nucleation, we used very large in- 
teraction ranges. We chose N — 8193 and R — 4096 
for the PDB model and N = 262144 and R ^ 13107 
for the PS model. These larger system sizes necessarily 
take a long time to run. For the PDB model, which even 
for smaller rangges has fairly long run times, this makes 
obtaining accurate statistics difficult. This difficulty is 
compounded by the fact that the spinodal is a critical- 
like point at which the correlation time blows up. This 
critical slowing down results in fairly large error ranges 
for X- 

We plotted x vs T simulation results in the insets 
of Fig. |3]d and Fig. [6|d. A divergeance appears as T 
approaches the spinodal temperature T^. This diver- 
gence appears power law when plotted on a log-log scale 
vs. Ts — T (Fig. [sjs and Fig. [6}d). Using a power-law fit, 
we calculated the susceptiblity exponents 7 « 0.7 ± 0.2 
and 7 « 0.50 ± 0.01 in the PDB and PS models, respec- 
tively. These exponents are consistent with our earlier 
theoretical results. 

C. Droplets and Growth Modes 

The critical droplet sits atop a saddle point hill in a 
free energy landscape. At this saddle point, the system 
has an equal probability of nucleating or returning to the 
metastable well. To determine the shape of the critical 
configuration, we first run our system until the strands 
have clearly separated. We then make 20 copies of the 
system and rerun them using the same initial conditions 
and random number sequence until some intervention 
time ti. At this time, each copy of the system is given a 
new sequence of random numbers that is different from 
each of the others. We then run the copies with the new 
random number sequences and calculate the percentage 
that still separate. The critical droplet is the configura- 
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FIG. 9; (a) Droplet profile and (b) initial growth in the PS FIG. 10: (a) Droplet profile and (b) initial growth in the PS 
model for TV = 128, i? = 1 and T = 365 K. The system uses model for iV = 4096, R = 205 and T = 470 K. The system 



a coarse-graining range Lcg ~ 5. 



uses a coarse-graining range Lcg = 205. 



tion that has approximately equal likelihood of growing 
into the stable phase or returning to the metastable well. 



In Figs. |9p, and 10 1, we plot the critical droplet profiles 
for the PS model with ranges R = 1 and R ~ 205, re- 
spectively. These profiles were taken at nucleation times 
obtained using the intervention procedure described pre- 
viously. For easier viewing, each profile has been coarse- 
grained by averaging over all bps within a length icGj 

i +Lc G 

(t>CG,i = X! 

i-LcG 

where (f) is replaced by a. The coarse-graining length Lcg 
is listed in each figure. Consistent with both our theo- 
retical results and results of previous studies [24142 7j . we 
obtain compact, large amplitude droplets for R = I and 
temperatures near the melting temperature, and diffuse, 
small amplitude droplets for R — 205 and temperatures 
near our theoretical prediction for the spinodal (note the 
change in scale). In addition to the droplet profiles, the 
growth modes were obtained by subtracting the critical 



profile from profiles at later times. The NN model shows 
an initial growth that is peaked at the surface of the 
droplet while the LR model give rise to growth modes 
that are peaked at the center of the droplet. These re- 
sults are consistent with the notion that nucleation in 
the lower temperature NN system occurs near a coexis- 
tence curve, while the LR system at a higher temperature 
nucleates close to a spinodal [24H27] . 

In Figs. [TT^ and [12^ , we plot the critical droplet pro- 
files for the PDB model with R — 1 and R = 205, respec- 
tively. As before, critical profiles were taken at nucleation 
times determined by intervention. The Morse well be- 
comes fiat in the range 1 < y < 2, so bps with separations 
greater than 2 can be considered open. The NN critical 
profile exhibits a compact open region with very large 
bp separations, i.e. yi > 10. Unlike the PS model, the 
region spans almost the entire length of the droplet (see 
discussion below). In contrast to the NN system, the LR 
PDB model produces diffuse critical droplets that differ 
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FIG. 11: (a) Droplet profile and (b) initial growth in the 
PDB model for iV = 128, J? = 1 and T = 380 K. 



little from the metastable background. These droplets 
are similar to those found using the LR PS model. Us- 
ing Eq. [25] with replaced by y and Log = 205, we 
coarse-grained the LR droplet profile to obtain a sharper 



image of its shape (red curve in Fig. 12 1. As with the 



PS model, we found growth modes for the PDB model by 
subtracting the critical profile from profiles at later times. 
The results are displayed in Figs. [TT] and [T2]d for vari- 
ous times after the critical droplet. Though noisy, one 
can easily see that growth occurs at the droplet 's edge 
in the NN system and at the droplet center in the LR 
system. These growth modes are consistent with those 
found in the PS models and those predicted by nucleation 
theory [24H22]- 

Classical nucleation theory predicts compact droplets 
greater than some critical size will grow from their edges 
to bring the system into the stable phase. This does not 
appear to be the case for NN PDB droplets. As noted 
earlier, the NN PDB droplet shown in Fig. [Tl] has the 
vast majority of 128 bps already open, with only about 
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FIG. 12: (a) Droplet profile and (b) initial growth in the PDB 
model for TV = 4096, R = 205 and T = 800 K. The droplet 
is shown in black. The red dashed curves are coarse-grained 
plots of the droplet and growth mode where the base pair 
separation has been averages over a coarse-graining length 
Lea = 205. 



20 bps still bonded. If the simulation is run using the 
same parameters but with different random noise, similar 
results are obtained, again with roughly only 20 bps re- 
maining bonded in the critical configuration. Moreover, 
if we increase the system size to iV = 256 but keep the 
temperature, range, and other parameters fixed, we find 
that the system again transitions to the unbound state 
with a critical configuration that has the vast majority 
of its bps unbound and roughly 20 still in the bonded 
state. Since the critical droplet grows with the size of 
the system, the NN model does not appear to be under- 
going classical nucleation. It is possible this is due to a 
finite size effect, but it is difficult to ascertain whether or 
not this is the case. In principle, we should be able to in- 
crease the system size until finite size effects go away, but 
in practice this is quite difficult. As mentioned earlier, 
the large flutuations that arise in the NN model make 
it difficult to see when metastable equilibrium has been 
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reached, so one must lower the temperature to produce 
longer metastable lifetimes. Increasing the system size 
makes these already very long runs even longer. Further- 
more, the range of interaction is much smaller than the 
system size, so it seems unlikely that this discrepency can 
be explained purely by finite size effects. 

If finite size effects do not cause the discrepency be- 
tween NN PDB droplets and those predicted by classical 
nucleation theory, then melting in the NN PDB model 
likely occurs via some mechanism other than nucleation. 
There is evidence to support this hypothesis. First, while 
the transition is sharp with a peaked in the heat capac- 
ity near the melting temperature, it has not yet been 
shown to be a true phase transition 3, 4 . This is in con- 
trast to the NN PS model, which includes the entropy of 
three-dimensional unbound loops to obtain a true phase 
transition when c > 2 |18) . Second, our simulations of the 
NN PDB model show recombination of the strands even 
after complete separation, suggesting again that melting 
is more likely a fluctuation than a first-order phase tran- 
sition. Furthermore, Van Hove argued that phase tran- 
sitions do not occur in one-dimension for systems with 
finite-range interactions |32j . Strictly speaking, his argu- 
ment does not apply here since it assumes the absence 
of external potentials. Indeed, one can easily see that 
one-dimensional phase transitions do occur by consider- 
ing a one-dimensional Ginzburg-Landau model with the 
(f)'^ polynomial expansion terms treated as an external 
potential. Unlike the Ginzburg-Landau model, the NN 
PDB does not have a second energetic well to represent 
the bound state. As shown in Section IV, the mean field 
PDB model does exhibit an entropic well caused by the 
nonlinear term in the coupling interaction, but the mean 
field model clearly contains LR interactions. It seems 
likely that the nonlinear coupling term is insufhcient to 
induce a phase transition via nucleation when the inter- 
action range is NN. In short, while the melting transition 
in the NN PDB model seems to be an activated process 
induced by a fluctuation, it does not appear to transi- 
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FIG. 13: Metastability for Heterogeneous Sequences. Plots 
are shown for the time evolution of (a) the = 128 PS model 
with R=l quenched to T = 365 K, (b) the iV = 128 PDB 
model with i? = 1 quenched to T = 380 K, (c) the iV = 4096 
PS model with 7? = 205 quenched to T = 465 K, and (d) the 
N = 4096 PDB model with R = 205 quenched to T = 805 K. 



tion via true nucleation. As one increases the range of 
interaction, the transition is sharper and can be more 
accurately described by nucleation. 

D. Heterogeneous Nucleation 

Biological DNA is heterogeneous, so if our results are 
to have biological significance they must hold for het- 
erogeneous systems. To test whether or not this is the 
case, we simulated the PS and PDB models for randomly 
generated sequences. Each bp was assigned a dissocia- 
tion energy of either Dj^t or Dqc with equal likelihood. 
As with our homogeneous results, we again monitor the 
time evolution of both models by plotting (y) and (<t) for 
R = 1 and R ~ 205 in Fig. [13] As before, each run is 
quenched from T — 300 K to a higher temperature whose 
value is chosen to give a reasonable run time. The value 
of the quenched T is given in the caption of each figure. 
As with homogeneous systems, we see a brief period of 
stability prior to spontaneous growth in both models for 
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FIG. 14: Coarse-grained dissociation energy, droplet profile, 
and initial growth in the PS model for A'^ = 128, R = 1 and 
T = 365 K. The system uses a coarse-graining range Lcg = 5. 



both ranges, indicating that the heterogeneities do not 
ehminatc metastabihty. 

For each of the runs in Fig. |13[ we determine the nu- 
cleation time using the intervention procedure described 
earher. For each model and range, we again plot the crit- 
ical droplet profile and growth modes (Figs. I4p7 ). In 
the top of each figure, we plot the coarse-grained dissoci- 
ation energy so that one can easily visualize how binding 
energy changes across the sequence. For each site i, the 
coarse-grained dissociation energy DcG,i is given by 



DcG, 



i+LcG 

E 

l-LcG 



D, 



(26) 



FIG. 15: Coarse-grained dissociation energy, droplet profile, 
and initial growth in the PS model for N = 4096, R = 205 
and T = 465 K. The system uses a coarse-graining range 
Lcg = 205. 



Though the results are nosier than the homogeneous case, 
we still observe that for both the PS and PDB NN models 
where quenches are restricted to be near the coexistence 
curve, melting is initiated by compact large-amplitude 
droplets. As before, the critical drolet of the NN PDB 
model spans most of the system. The growth modes 
of both NN models feature peaks at the surface of the 
droplet. In both models, LR interactions allow the sys- 
tem to reach metastabihty at higher temperatures. This 
brings the system closer to an apparent pseudospinodal, 
which gives rise to diffuse critical droplets whose ampli- 
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FIG. 16: Coarse-grained dissociation energy, droplet profile, 
and initial growth in the PS model for A'^ = 128, R = 1 and 
T = 380 K. 



tude is close to the metastable background. Both LR 
systems grow from the center of the droplet. These re- 
sults suggest that the inclusion of random heterogeneities 
does not have a major effect on the qualitative features 
of the transition. 

The formation and growth of a critical droplet will dif- 
fer from seqeunce to sequence and even from run to run 
for a particular sequence of bps. Still, it is worthwhile to 
analyze the particular pathways that the heterogenous 
systems reported here have taken. First, one will note 
that droplets of both models tend to be peaked around 
areas where the binding strength is weak. For example. 



FIG. 17: Coarse-grained dissociation energy, droplet profile, 
and initial growth in the PS model for N = 4096, R = 205 and 
T = 850 K. The red dashed curves are coarse-grained plots of 
the droplet and growth mode where the base pair separation 
has been averages over a coarse-graining length Lcg ~ 205.. 



the NN PS model features a dip in the binding energy 
around bp index 70, which also corresponds to the peak 
of the droplet. Similarly, the critical droplet in NN PS 
model features two peaks around bp indexes 50 and 100, 
which correspond to dips in the binding energy. Partic- 
ularly noticible is the dip in the critical droplet near bp 
index 75 that appears between these two peaks. This dip 
corresponds to a stronger than average binding strength. 
A similar dip in the LR PDB droplet occurs between bp 
indexes 1000 and 2000, at which there is a maximum 
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in the binding energy of the strands. While the criti- 



cial droplet of the LR PS model shown in Fig. 15 does 
not correspond to a large dip in the binding energy, the 
growth mode appears to develop a shoulder on the left 
side of the droplet at which there is a prominent dip 
in the binding energy. These results suggest that while 
modest amounts of inheterogeneities do not completely 
distort the qualitative character of the transition, they 
do make melting pathways at weaker sites more probable 
than pathways through strongly bound bps. This would 
imply that heterogeneities cause only minor changes in 
the droplet structure but cause appreciable changes in 
the nucleation rate. 

VI. SUMMARY AND DISCUSSION 

We have extended the PS and PDB models of DNA 
to include variable-range interactions and have observed 
nucleation in the both the original and long-range mod- 
els. In the PS model, we constructed a ^''-field theory by 
enumerating all possible spin configurations of a coarse 
grained region, calculating the average entropy of the re- 
gion due to both loops and coarsening, and assuming 
that this average was equivalent to the total entropic 
contribution from the region. While this field theory is 
not expected to produce accurate results near the coexis- 
tence curve where the coarse-graining size is smaller than 
the largest loops, it is useful in calculating droplet pro- 
files and growth modes at higher temperatures where the 
droplets become more diffuse. In addition, the field the- 
ory predicts a divergent susceptibility near the spinodal 
with a pseudocritical exponent 7 0.5. 

We constructed a self-consistent mean field theory of 
the PDB model from which we calculated x{T) and ob- 
served a divergence as T characterized by the spin- 
odal exponent 7 0.5. 

Simulations of both models for a variety of ranges 
show some signature characteristics of nucleation includ- 
ing metastability prior to growth and a constant tran- 



sition rate. Critical droplets and growth modes were 
measured both near the coexistence curve and near the 
spinodal in systems with short and long range interac- 
tions, respectively. We find that the NN PS model pro- 
duces compact large amplitude droplets similar to those 
of classical nucleation theory, while the NN PDB model 
produces droplets that span nearly the entire system size. 
Both NN models appear to grow from the droplet surface. 
Both the PDB and PS models exhibit diffuse small ampli- 
tude droplets when given LR interactions and quenched 
near the theoretically-predicted spinodal. As expected, 
we find compact classical droplets growing predominantly 
at their surface and diffuse spinodal droplets growing 
from their center. The presence of heterogeneous se- 
quence does not appear to appreciably alter the above 
results. By measuring the fluctuations, we calculated 
x(T) and observed a divergence at large T characterized 
by the pseudo critical exponent consistent with our the- 
oretical results. 

The consistency of the critical exponents in both mod- 
els suggests that they may be in the same universality 
class for spinodals. By measuring the susceptibility near 
the melting temperature, it may be possible to experi- 
mentally determine how accurate a mean field depiction 
of DNA is and whether or not long-range interactions 
play a significant role in determining the qualitative char- 
acter of nucleation observed. In particular, the presence 
or absence of long-range interactions may determine the 
most effective way for biological enzymes to mechani- 
cally denature DNA locally, and an understanding of the 
qualtitative difference between the long and short range 
systems may be valuable to researchers experimenting on 
these enzymes. 
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